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Manipulating fluids at the nanoscale within networks of channels or chemical lanes is a crucial 
challenge in developing small scale devices to be used in microreactors or chemical sensors. In this 
context, ultra-thin (i.e., monolayer) films, experimentally observed in spreading of nano-droplets or 
upon extraction from reservoirs in capillary rise geometries, represent an extreme limit which is of 
physical and technological relevance since the dynamics is governed solely by capillary forces. In this 
work we use kinetic Monte Carlo (KMC) simulations to analyze in detail a simple, but realistic model 
proposed by Burlatsky et al. A, % f° r the two-dimensional spreading on homogeneous substrates 
of a fluid monolayer which is extracted from a reservoir. Our simulations confirm the previously 
predicted time-dependence of the spreading, X(t — > oo) = Ayft, with X(t) as the average position 
of the advancing edge at time t, and they reveal a non-trivial dependence of the prefactor A on 
the strength Uo of inter-particle attraction and on the fluid density Co at the reservoir as well as 
an [To-dependent spatial structure of the density profile of the monolayer. The asymptotic density 
profile at long time and large spatial scale is carefully analyzed within the continuum limit. We 
show that including the effect of correlations in an effective manner into the standard mean-field 
description leads to predictions both for the value of the threshold interaction above which phase 
segregation occurs and for the density profiles in excellent agreement with KMC simulations results. 

PACS numbers: 68.15.+e, 68.37.-d, 81.15.Aa 



I. INTRODUCTION 

In the context of microfluidics, wetting phenomena at 
the micro- and nano- meter scale 0, SIS 00] are relevant 
for applications like microreactors or chemical sensors, 
for which a crucial challenge is the transport of liquid 
to networks of channels or chemical lanes, as well as its 
precise manipulation within such a network 0, H Il0| . 
Since at this small scale the liquid-substrate interaction 
is important, the flow of thin films may be eventually 
controlled by engineering the physico-chemical properties 
of the substrate, thus opening the road for applications 
which do not have an equivalent at the macroscopic scale 

mm us. 

Although the existence of very thin precursor films has 
been long ago evidenced by the studies of Hardy 
only recent experiments on liquid spreading on atomi- 
cally smooth surfaces EM EH E3, EM HH H^l , 
performed with volumes of the order of nano-liters, have 
clearly shown by means of dynamic ellipsometry or X- 
ray reflectivity measurements that one or few precursor 
films with molecular thickness and macroscopic extent 
advance in front of the macroscopic liquid wedge of the 
spreading drop. The liquids used were low-molecular- 
mass polymer oils which behave as non-volatile liquids, 
and experiments performed both for spreading of nano- 
droplets and for capillary rise geometries have established 
that the linear extent X(t) of the precursor film grows in 
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time as X(t — ► oo) ~ At a . The exponent a = 1/2 seems 
to be independent of the nature of the liquid and of the 
substrate, of the geometry, and of the volume of droplet 
as long as the droplet is not emptied and constitutes a 
reservoir for the extracting film; only the prefactor A de- 
pends on these parameters. 

Several theoretical models have been proposed (see 
Refs. S S El El and references therein) and 
an impressive number of Molecular Dynamics (MD) and 
Monte Carlo (MC]numerical simulations have been per- 
formed (see Refs. (SH^miliEmiillHlil and ref- 
erences therein) in order to understand the mechanisms 
behind the extraction of precursor films and to explain 
the time dependence of the spreading. 

The hydrodynamic model of de Gennes and Cazabat 
[23| assumes a layered structure of the droplet, each layer 
being a two-dimensional incompressible fluid, in which 
vertical transport is possible only at the edges of the lay- 
ers. The model leads to the correct time dependence for 
the advancing layers, but, as pointed out in Refs. 0, [24| . 
and |34j , it is debatable if this hydrodynamic description 
holds at the molecular level and can be directly applied 
to ultra-thin films. 

A different approach, along the line of earlier work on 
activated kinetics by Cherry, Holmes, Blake, and Haynes 
[35l l3o| , consists of a microscopic description for the thin 
liquid films in terms of lattice gas models for interacting 
particles. One such model is the two-dimensional driven 
Ising model recently proposed by Abraham et al. |26j| . 
Using kinetic Monte Carlo (KMC) simulations it was 
shown that in this model the transport of mass occurs via 
a second layer, and a particle-hole diffusion equation was 
used to show that the model leads to correct predictions 
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(confirmed also by simulations) for the time dependence 
of spreading. The model predicts a uniform density and 
a compact first monolayer, in close resemblance of the 
incompressible layers of the hydrodynamical model men- 
tioned before. 

For the case that the precursor consists of a single 
monolayer, a lattice gas model of interacting particles has 
been proposed by Burlatsky et al. This model, which 
allows mass transport from the reservoir to the advanc- 
ing edge only inside the monolayer, has been extended to 
the more general situation of relaxation of a monolayer 
initially occupying a half-plane without a reservoir by 
Oshanin et al. 2]. Based on several mean- field assump- 
tions, among which the strongest is the replacement of 
the inter-particles attraction in the "bulk" by an "effec- 
tive force" acting on the advancing edge, the authors have 
been able to derive the t 1 / 2 -law for the spreading and to 
calculate the dependence of the prefactor on the fluid- 
fluid interaction parameters. In contrast to the other two 
models, in this case the density in the monolayer depends 
significantly on the distance from the reservoir. Although 
it is reasonable to expect that neglecting the attractive 
particle-particle interactions in the "bulk" should not af- 
fect the time dependence ~ t 1 / 2 , one can expect that 
the behavior in the presence of attractive interactions is 
much richer (see, for example, the recent KMC simula- 
tions results of Lacasta et al. [31} for a closely related 
one-dimensional model) . 

All three models mentioned above recover the correct 
time-dependence of spreading, but it is unlikely that one 
can discriminate between them via experimental tests 
based on their predictions for the corresponding pref- 
actors because these include in each case a number of 
parameters whose connections with experimental quan- 
tities are not clear. However, we have already pointed out 
that these models lead to qualitatively different predic- 
tions with respect to the shape of the emerging density 
profiles (constant in the first two cases, spatially vary- 
ing in the last) which are, in principle, experimentally 
accessible. 

In this work we analyze in detail the density profiles of 
this last model. We present results of KMC simulations 
on a square lattice of a model for spreading of a liquid 
monolayer closely related to the model in Refs. and Q. 
Our choice for KMC simulations of a lattice gas model 
is motivated by the fact that we are interested in the 
asymptotic (large spatial and temporal scales) behavior, 
a regime which as yet cannot be explored using Molecu- 
lar Dynamics simulations because of extensive computing 
resources needed to simulate spatially large systems and 
unreasonable large CPU times required to simulate real 
times even in the order of /z-seconds. In contrast to the 
previous work mentioned above, we shall explicitely con- 
sider the asymmetry of the jump rates in the bulk, at the 
expense of being able to measure the prefactor A from 
the simulations but not to predict it analytically. Our 
results show a non-trivial dependence of the prefactor A 
on the strength Uq of the inter-particle attraction and on 



the density Co at the reservoir. The asymptotic spread- 
ing behavior at long time and large spatial scale of the 
transversally averaged density profile is analyzed within 
a continuum limit. We show that the model predicts 
qualitatively different structures for the experimentally 
accessible density profiles along the spreading direction 
above and below a threshold value for the ratio between 
the fluid-fluid interaction and the thermal energy. In- 
cluding the effect of correlations in an effective manner 
into the standard mean-field description, we find excel- 
lent agreement between the theoretical predictions and 
the KMC results. We conclude the paper with a sum- 
mary and discussion of the results. 



II. DEFINITION AND DISCUSSION OF THE 
MODEL 

As mentioned in the Introduction, a simple microscopic 
model for the dynamics of a fluid monolayer in contact 
with a reservoir was proposed in Refs. p| and 0- Al- 
though we use here this lattice gas model of interacting 
particles with only slight modifications, for clarity and 
further reference we describe, motivate, and comment on 
the defining rules as follows. 

(a) The spreading geometry is rectangular (x — y plane) 
and the substrate is homogeneous. The half-plane x < 
is occupied by a reservoir of particles (three-dimensional 
bulk liquid) at fixed chemical potential which maintains 
at its contact line with the substrate, positioned at the 
line x = 0, an average density Co (defined as number 
of particles per unit length in the transversal y direc- 
tion). For the case of capillary rise, the reservoir would 
correspond to the liquid bath and the line x = to the 
edge of the macroscopic meniscus. It is assumed that 
the only role of the reservoir is to maintain Cq constant, 
and thus to feed the monolayer which is extracted, but 
there is no flow of particles from the reservoir to "push" 
the film. The parameter Co is expected to be related 
to the difference in the free energy per particle A_F be- 
tween a fluid particle in bulk liquid, i.e., inside the three- 
dimensional reservoir, and one on the surface of the sub- 
strate at x = 0, and thus it is a measure of the wettabil- 
ity of the substrate by the liquid |2l|. A general, explicit 
form for the relation between Cn and AF is not available, 
but for a qualitative picture 0, an argument based 
on Langmuir-type adsorption may be used to estimate 
Co ~ aC reservo i r [l — exp(— /3AF)], where a is the area 
per adsorption site and C reservo i r is the density (number 
of particles per unit volume) in the reservoir. At time 
t = 0, the half- plane x > is empty. 

(b) The substrate-fluid interaction is modeled as a pe- 
riodic potential forming a lattice of potential wells with 
coordination number z and lattice constant a. The parti- 
cle motion proceeds via activated jumps between nearest- 
neighbor wells; evaporation from the substrate is not al- 
lowed. We assume that the dynamics of the activated 
jumps can be described by the classical reaction-rate the- 
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ory, i.e., the activation barrier Ua is significantly larger 
than the thermal energy kgT, where kg is the Boltz- 
mann constant and T the temperature, and the coupling 
to the substrate is large enough such that in crossing 
the barrier all the kinetic energy of the particle is dis- 
sipated [23. The barrier Ua determines the jumping 
rate Vt = vq exp[— UA/ksT], where vq is an attempt fre- 
quency defining the time unit. We note that for a two- 
dimensional homogeneous, isotropic substrate and a reg- 
ular (square, triangular, honeycomb, etc.) lattice struc- 
ture, this jumping rate can be absorbed into the one- 
particle diffusion coefficient Dq = Via 2 /A on the bare 
substrate |3j|. Therefore, in this case Ua is an irrele- 
vant parameter in the sense that it can be incorporated 
either into the choice of the unit of time as f2 _1 or into 
that of the unit of length as y/ Dq/i/q. For the rest of this 
work we consider a square lattice (z = 4); a qualitatively 
significant dependence of the results on the lattice type is 
not expected. This expectation is based on correspond- 
ing results obtained from test simulations on a triangular 
lattice (z = 3). 

(c) The pair interaction between fluid particles at dis- 
tance r is taken to be hard-core repulsive at short range, 
preventing double occupancy of the wells, and attractive 
at long range, —Uo/r 6 for r > 1, resembling a Lennard- 
Jones type interaction potential. Here and in the follow- 
ing r is measured in units of a so that Uq denotes the 
strength of the interaction energy. The absence of dou- 
ble occupancy leads to an a priori removal of thickening 
of the film as a possible relaxation mechanism, which is 
not meant to imply that we consider it irrelevant. We 
have decided to disregard this mechanism here since it 
would have significantly increased the complexity of the 
problem. Thus we leave the issue of film thickening open 
for further research. 

(d) As we have mentioned in (a), the motion proceeds 
via activated jumps between nearest-neighbor wells, the 
activation barrier for any jump being Ua- The selec- 
tion of the nearest-neighbor well, i.e., the probability 
p(r — > r'; t) that a jump from location r will be directed 
toward the location r', is biased by the fluid-fluid energy 
landscape and is given by 



for the numerical simulations, 
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exp{f[fr(r;f)-fr(r';t)]} 
Z(r:t) 
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where Z(r;t) = ^ expj -[U(r;t) - U(r';t)} \ is 
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the normalization constant and 1//3 = fcsT, 
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and rj(r'\ t) 6 {0, 1} is the occupation number of the well 
at r' at the time t. We note that, after canceling the 
common factor exp[[3U(r; t)/2], the expression Eq. QJ 
may be re-written in a form which is somewhat simpler 



exp 



p(r -> r';t) 



~U{r';t) 



E exp 

r',|r'— r\ — 1 



(3) 



the dependence on r being retained because the summa- 
tion is carried out over the neighboring locations of r. 
We also note that the summation in Eq. © has been 
restricted to three lattice units for computational conve- 
nience. This corresponds to the cut-off generally used in 
Molecular Dynamics simulations for algebraically decay- 
ing Lennard-Jones pair-potentials. 

The expression Eq. for the probability that a cer- 
tain direction is chosen for jumping deserves further dis- 
cussion. For a particle located at r, it follows from the 
definition of p(r — > r';t) (Eq. QJ) that the rates 



= Vlp(r -> r';t) 



(4) 



for the transitions from r to neighboring locations r' sat- 
isfy 



r' t*| = 1 



(5) 



Thus the total rate of leaving a location for any given 
particle at any given location is determined only by the 
fluid-solid interaction characterized by Ua, it is time- 
independent, and it equals Q. Therefore, the fluid-fluid 
interaction will influence only the choice of a direction for 
the jump, but not the jumping frequency, and in the dy- 
namics there will be only one relevant microscopic time 
scale, f2 _1 , which is dictated by the solid-fluid coupling. 

The choice p(r -> r';t) cx exp{§ [U(r; t) - U(r';t)}} 
is motivated by the following. If one disregards the 
reservoir and considers a system with a given volume 
and a given number of particles, the change in the 
total fluid-fluid energy (no double occupancy) Ut — 
l/2^2 r rj(r;t)U(r;t) due to a change in configuration 

{n{r) = I, r,(r') - 0} -> { V (r) = 0, V (r') = 1}, (6) 

where \r — r'\ = 1, is given by AUt = U(r';t) — U(r;t), 
with U(r'; t) calculated for the final configuration and 
U(r;t) for the initial one. Then a simple choice of tran- 
sition rates to' r _ >r ,. t which satisfies detailed balance with 
respect to the equilibrium distribution 



where Z 



E 

all {r](r;t)} 



V = exp(-0U t )/Z, 
exp(-f3U t ), is 



(7) 



u>' r ^ r ,. t =nexp{-pAU t /2). (8) 
If the transition rates would be chosen according to the 
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expression Eq. JSJ above, then the fluid-fluid interaction 
would effectively change the activation barrier and lead 
to a whole spectrum of microscopic time scales. Normal- 
izing these local transition rates (Eq. JSJ) by the (local) 

total transition rate id r ^ r ,. u a decoupling re- 

r',|r'— r| — 1 

suits between an activated dynamics determined by the 
solid-fluid interaction and a weak perturbation due to the 
fluid-fluid interaction, and the choice given in Eq. Q for 
the bias probability p(r — > r'; t) is obtained. One should 
note that this decoupling is obtained at the expense that 
the transition rates uj r ^ r i- A defined by Eq. (0J do not 
satisfy detailed balance with respect to the distribution 
in Eq. (0, although the deviations, which are equal to 
Z(r;t)/Z(r';t) due to 



Z(r;t) 
Z(r';t) 



exp(-(3AU t ), 



(9) 



are expected to be very small. From this point of view, 
the dynamics defined by Eq. is that of an asymmetric 
exclusion process |40|. but with a position- and time- 
dependent bias. 

(e) As defined by the rules (a)-(d), the model cor- 
responds to mass transport from the reservoir into a 
two-dimensional vacuum so that a phase with very low- 
density, due to two-dimensional evaporation, will form in 
front of the advancing monolayer. The emergence of this 
low-density phase poses problems in that its long time 
dynamics, which is of ideal gas type, mixes with that 
of the following-up "compact" film and leads to serious 
difficulties in defining the advancing edge of the mono- 
layer. This problem has been encountered earlier also in 
three-dimensional simulations IH HI |H HI E3 and, 
in general, it has been overcome by replacing the sim- 
ple particles by connected chains mimicking polymers. 
Although this approach is straightforward it is not very 
appealing, neither from a theoretical point of view (an 
analytical approach becomes at least cumbersome, if not 
intractable), nor from a computational one (the memory 
and CPU requirements for sufficiently long simulation 
runs for large enough systems are unreasonably high). 

Therefore, we have adopted a different approach. We 
define the advancing edge T t of a monolayer configuration 
at time t as the set of the most advanced particles in each 
line y = const for this configuration (see Fig.^l. We elim- 
inate two-dimensional "evaporation" by imposing the fol- 
lowing additional constraint: moves from sites r G T t to- 
ward sites r' ahead of T t for which \U(r';t)\ < U c , where 
U c > is a fixed threshold value, are rejected. This corre- 
sponds to requiring a given minimum number of particles 
in the neighborhood \r\ < r c of any of the components 
of T t . The results presented in this paper correspond to 
simulations with U c = Uq/3 6 , i.e., r c = 3, in other words 
to the requirement that in the disk \r" — r'\ < 3 centered 
at r' there is at least one more particle in addition to the 
one attempting the jump r — > r' (see also, c.f., Fig.[5j|. 

The above constraint is close in spirit to the "ef- 
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FIG. 1: Schematic drawing of a typical configuration {r/(r; t)} 
of a monolayer spreading on a rectangular lattice (viewed un- 
der an oblique angle). The particles denoted as open circles 
occupy the edge of the reservoir x = with Co(t). Black cir- 
cles denote the particles in the bulk of the monolayer whereas 
gray circles denote particles at the advancing edge Ft- Also 
indicated is the average position X{t) of the advancing edge. 
Here Co(t) and X(t) correspond to averages over y for a given 
realization and not to averages over runs (for which we use the 
same notation in the main text). There are periodic boundary 
conditions in the y— direction. 

fective boundary-tension" idea used in Ref. in 
which the attractive interactions have been neglected 
except for particles on the advancing edge for which 
a constant asymmetry in the jumping rates "away" 
and "toward" the reservoir was imposed. Rule (e) 
provides a simple and convenient way of controlling 
the rate of two-dimensional evaporation. For example, 
setting U c = corresponds to fully unconstrained 
dynamics, while replacing the rejection procedure with 
an "acceptance rate" will allow for a continuous tuning 
of the evaporation rate through the acceptance rate. 
We note that, physically, the model defined by the rules 
(a)-(e) could be used to study also expansion into 
an already present vapor phase instead of expansion 
into vacuum. In the presence of a vapor phase there 
would be an average occupancy of the sites in front 
of r t , and thus some of the jumps from T t would 
be rejected due to the hard-core repulsion, which is 
an effect similar to an acceptance rate as discussed above. 



III. KINETIC MONTE CARLO SIMULATIONS 

We have carried out KMC simulations of the model 
defined in Sect. II using square lattices with widths L y 
of 200 or 500 lattice units, periodic boundary conditions 
along the transversal (y) direction (appropriate for simu- 
lating an infinitely wide substrate) , and an activation en- 
ergy [3Ua — 3.5. Some simulation runs have been carried 
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out using lattices with smaller widths in order to check 
finite-size effects. We have found that for widths larger 
than 100 lattice units there is no detectable influence of 
the width value on the quantities we have measured in 
these simulations. The length of the lattice in the x direc- 
tion has been chosen to be L x = 1000 lattice units, with 
the possibility of changing it dynamically in the course 
of the simulation if necessary, i.e., if T t intersects the 
line x = 1000; however, this situation was not encoun- 
tered in any of the simulations we have carried out. We 
note here that in the experiments mentioned in Sect. I 
El 13 IH 13 13 II El typical values for the diffusion 
coefficient were estimated from the spreading rate to be 
of the order of 10 -11 — 10 -9 m 2 /s. If we take the lattice 
spacing as a ~ 10 nm, i.e., of the order of the lateral size 
of a PDMS coil |22(, then the above values for the dif- 
fusion coefficient imply typical values for the frequency 
n of the order of 10 5 - 10 7 s" 1 , and thus for (3U A = 3.5 
typical values for v of the order of 10 6 — 10 s s _1 . 

For the simulations we have used a variable step con- 
tinuous time kinetic Monte Carlo algorithm 0, OH 
which is described in Appendix A. One step in the Monte 
Carlo simulation proceeds as follows. At time t, a par- 
ticle from the film (x > 0) is selected at random. The 
time is incremented with At (the time at which a jump 
attempt with sufficient energy for leaving the well will 
occur), where At is a random variable distributed ac- 
cording to P( At) = Nil exp and 
N is the number of particles in the film at time t (so that 

(At)p = j^q)- The direction for the jump is chosen at 

random with probabilities weighted according to Eq. • 
If the destination site is empty, the jump takes place; if 
not, the jump is rejected. Exchange between the reser- 
voir [x < 0) and film (x > 0) is subject to the additional 

constraint that the density Co(t) = — \^T](x = 0,y;t) 

L V y=l 

on the line x = fluctuates narrowly around a given 
value Cq and proceeds in the following manner. Moves 
from a; = to i = -1 are allowed if Co{t) is maintained 
within the interval [Cn(l — S), Co(l + S)], where the am- 
plitude 8 has been fixed to 10~ 2 . If this condition is 
satisfied, the particle is considered to become part of the 
reservoir and is removed; if not, the move is rejected. In 
case of moves from x = to x = 1, if the density on 
x = would decrease below Co(l — 6), then after the 
move a new particle is added on an empty site (chosen 
at random) on the line x = 0. Similarly, in case of moves 
from i = 1 to i = a particle is removed (at random) 
from the line x = Q if the density would increase above 
Co(l + S). The time is not incremented upon adding or 
removing particles, corresponding to the reasonable as- 
sumption that the equilibration of the reservoir is very 
fast. In order to compute the potential energy at a desti- 
nation site on the line x = — 1, needed for the evaluation 
of the weight probabilities for jumping of a particle on the 
line x = 0, in the beginning of the simulation particles 



are placed at random on the lines —4 < x < — 1 such that 
the average density on these lines is Co, and this configu- 
ration is kept unchanged during the simulation run. Due 
to these procedures one does not have to consider the 
dynamics on the lines x < —1, i.e., in the reservoir. In 
order to have sufficient fluctuations in Co(t) on the line 
x = 0, fluctuations which mimic the stochastic nature of 
the exchange of particles between the reservoir and the 
film, the width L y has to be large enough such that the 
amplitude <5 translates into a reasonable number of sites. 
This is the reason why we have used a large value for the 
width; for example, at Co = 0.8 and L y — 500, S = 10~ 2 
translates into 4 sites. 

All the "measured" quantities have been averaged over 
a number of independent simulation runs ranging from 
10 to 50, a value of 50 runs being used in most of the 
cases. These runs differ from each other both with re- 
spect to the initial configuration {i](x — 0,y;t = 0)} 
and the subsequent sequence of jumps. The observables 
of interest are defined below. The density p(r;t) is de- 
fined as p(r;t) — (r](r;t)), where (•••) means average 
over different KMC runs. Due to the symmetry of the 
model, the density profile p(r; t) in the limit of infinitely 
many runs is independent of y, while in the average over 
a finite number of runs random, uncorrelated fluctua- 
tions (whose amplitude decrease with increasing number 
of runs) occur along the y direction. These fluctuations 
are suppressed by measuring the transversally averaged 

/ 1 Lv \ 
density C(x, t) = ( — 77(2;, y:t) \ , and thus it is ex- 

\ v y=i I 
pected that C(x,t) ~ p(r;t), with strict equality for in- 
finitely many runs. The average position of the advancing 



edge of the monolayer is defined as X(t) — 




For the case U c = Uo/3 e (as used for the actual sim- 
ulations), two-dimensional evaporation is negligible and 
X(t) (which we shall also call front line) is a good mea- 
sure for the actual advancing edge of the monolayer. We 
note here that in the case when no constraint is imposed 
to prevent two-dimensional evaporation (U c = 0), the 
front line may be defined as X(t) — (x(t)) |44j . where 
x(t) is the most advanced line corresponding to a given 

1 L 

(smallest measurable) density C = — r)(x, y; t). Al- 

L v y=l 

ternatively, one may follow the time dependence of the 
mass of the film |37j , or use a percolation type definition 
for the precursor as the set composed of all the parti- 
cles which are connected along nearest-neighbor bonds 
with the reservoir, the boundary of this cluster defining 

r t M- 

Snapshots of typical density profiles during spreading 
are shown in Fig. [3 for the cases (a) Wq — 1.4 and (b) 
Wo = 0.8, where we have introduced the notation 

Wo - PU . (10) 
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FIG. 2: Typical density profiles p(x,y;t) (obtained by averaging over 50 KMC runs) for Wo = (3Uo = 1.4 (row (a)) and 
Wo = 0.8 (row (b)) at times (left to right) t = 2 x 10 5 , t = 10 6 , and t = 2 x 10 6 for Co = 1.0 and L y = 500. Time is measured in 
units of Ug = VC 1 exp(— /3Ua) with [3Ua = 3.5, distances are measured in lattice units, and spreading occurs in a;— direction. 
The color coding (shown on the right) is a linear function of density. 



These density profiles reveal already the qualitative dy- 
namical behavior. It can be seen that, as expected, the 
monolayer is homogeneous in the y (transversal) direc- 
tion, while along the x (spreading) direction there are 
significant density variations. As intuitively expected, 
the spreading of the edge of the monolayer is faster for 
smaller Wq, i.e., higher temperature (at a given interac- 
tion strength Uq) or smaller inter-particle attraction (at 
a given temperature T). In addition, one observes qual- 
itatively different dynamics as revealed by the abrupt 
change from high to low density for the large value of 
Wq compared to the smooth and broad decrease for the 
small value of Wq. This is accompanied by a different 
dynamics of the regions with moderate to low density in 
the two cases. While for large value of Wq (panel (a)) the 
range of densities p ~ 0.5 (the green band) and that of 
low densities p < 0.2 (the light blue band) maintain rela- 
tively constant widths and advance with similar rates, for 
smaller Wq (panel (b)) the regions become broader with 
increasing time and clearly the low density (light blue) 
has a larger rate of advancing. These features will be 
discussed more in the following sections, where a quanti- 
tative characterization of the spreading as a function of 
the parameters Wq and Cq is presented. 



IV. TIME-DEPENDENCE OF SPREADING 
AND ASYMPTOTIC SCALING OF DENSITIES 

A. Time dependence of spreading 

We start our analysis of the dynamics of spreading 
by studying the time-dependence of the position X(t) 
of the front line. As noted in the Introduction, we ex- 
pect that the asymmetry in the jumping rates due to 
the inter-particle interactions (Eq. QJ) will not affect the 
previously (without inter-particle interactions) predicted 
X(t) r-j \ft time dependence as long as the attractive 
interaction is not too strong. Figure shows the time 
dependence of the scaled front line, X(t)/s/Dot, for sev- 
eral values of the interaction parameter Wq and of the 
reservoir density Cq. 

One can see that at low and intermediate values of 
the interaction strength (Fig. 01 a)) the y/t behavior is 
recovered independent of the value Cq of the density in 
the reservoir. However, for strong attractive fluid-fluid 
interaction (Fig. Etb)) and low densities Co, the time 
dependence of X(t) clearly deviates from the y/t behav- 
ior, the latter being obtained only for high densities Cq. 
Since the decreasing trend shown by the data in Fig. |3{b) 
may be due to either spreading at a rate slower than 
or to the fact that there is no extraction of a macro- 
scopic film from the reservoir, we have also analyzed the 
time dependence of the total mass of the film extracted, 



7 



2.0 



0.0 



10 



ii 
o 
(hi 



^□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□c 

^qqOOCKXXXXXXXXXXXXXXXXXXXXXXXXXXXXiOOO 

" > <><> ooooooo«ooooooooooooooooooooooo©ooo 
^W =0.4 .»W = 1.O 



— — — — 1 



5~ "♦♦»*♦♦♦»♦♦♦♦♦♦♦♦♦♦»»>♦*«»*»**>>♦♦»♦♦*♦ 



Co = 1.0 : n, ■ 

Co = 0.8 : O, • 

C = 0.6 : 0, ♦ 



(a) 



time t 



2x10° 



10' 



W„ = 1.6 



(M C = 1.0:D 
v ' C = 0.8 : O 
C = 0.6 : 

□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□ 



10 H 



10' 



□ o o o o o 
o o o o o 



10 



time t 



10 




4& 



time t 



2x1 7 



FIG. 3: Front line position (divided by \/Dot) as a function 
of time for (a) Wo = 0.4 (open symbols) and Wo = 1.0 (filled 
symbols), and (b) Wo = 1.6, for reservoir densities Co = 1.0 
(squares), Co = 0.8 (circles), and Co = 0.6 (diamonds). The 
inset (on logarithmic scales) compares the time dependence 
of the (unsealed) mass M(t) of the film for Wo = 1.6 and for 
the same reservoir densities (symbols) with the t 1 ^ 2 behavior 
(dashed line). In (b) for Co = 0.8 and Co = 0.6 the scaled 
front line position decays ~ t — 1 / a for large t. Here and in the 
following time t is measured in units of v^ 1 = exp(— (5Ua) 



with (3Ua = 3.5, so that 



exp(-yU A ). 



M(t) = ( ^ )• As can be seen in the inset in 

\x>0,y I 

Fig. Efb), for large Co the time dependence of M{t) is 
very well described by \fi, as expected while for low 
values Co the mass M(t) shows a clear saturation to a 
small constant value (which corresponds roughly to a po- 
sition of the front of ca. 10 lattice units), thus indicating 
that actually there is no macroscopic film extracted from 
the reservoir. At early times fluctuations lead to the 
extraction or leakage (with a very fast decreasing rate 
of extraction) of a small number of particles from the 
reservoir, but the spatial extension X{t) of the film in 
this case remains microscopicly small and becomes time 
independent for f > 1. Therefore, as a function of the 
interaction strength Wo there is a transition from a "sub- 



strate covering" state at low values Wo, in the sense of 
extraction of a film with macroscopic lateral extension 
in the spreading direction independently of the density 
value Cq in the reservoir edge, i.e., a film which spreads 
according to X(t oo) ~ \fi, to a "non-covering" state 
at large values Wo, in the sense that a macroscopic film is 
extracted only for sufficiently large densities Co (eventu- 
ally for none if Wo is sufficiently large) . Results similar to 
the ones shown in Fig. El from simulations performed for 
a broad range of parameters values (0.4 < W < 1.6 and 

0.1 < Co < 1.0) indicate that the values W (cmj) (C ) for 
which this change in behavior occurs are bounded from 

below by 1.0 < Wq COV \Co) ■ We shall return to this point 
in the next but one paragraph and during the discussion 
of the continuum limit. 

The results presented in Fig. [3] also show that in case 
of spreading the time independent dimensionless prefac- 
tor A in X(t — > oo) = A^/Dot depends on both Wo 
and Co- From the curves X(t)j \fD~o~i one can estimate 
A = lim X(t)l \J Dot by fitting the data in the range 

t — >oo 

t ^ 1 (in practice the data in the last tenth of the 
time interval available) with a constant. The results 
for A(Co,Wo) are shown in Fig. Ufa). Here we use the 
convention that in the case where there is no macro- 
scopic film spreading in the sense explained above, i.e., 
X(t)/y/D~oi decreases in time and M(t) has a time de- 
pendence t ai with a.\ significantly smaller than 1/2 (in 
practice ct\ < 0.4), the value of the prefactor is assigned 
to be A = 0. As expected, A{Cq,Wq) is an increasing 
function of Co for fixed Wo and a decreasing function 
of Wo for fixed Co- However, the functional dependence 
is not simple, and one can easily notice a change in the 
shape of A(Cq,Wo) for Wq close to the value 1.0. For 
values Wo > 1, the curve shows a plateau over a range 
of densities Co which increases with increasing Wo, while 
for values Wo < 1 the prefactor A is a strictly increas- 
ing function of Co- This property is of a different type 
than the transition from covering to non-covering dis- 
cussed above, because it involves a change in the depen- 
dence of the prefactor A on the density Co while the 
spreading law X(t) ~ t 1 / 2 holds. This change in be- 
havior emerges (as we shall discuss in more detail in the 
next section) as a consequence of the competition be- 
tween the diffusive motion driven by the concentration 
gradient and the clustering tendency (opposing the con- 
centration gradients) driven by the inter-particle attrac- 
tion. This competition leads to instabilities at certain 
density values if the attractive interaction is sufficiently 
strong. From Fig. Ufa), the value of the threshold inter- 
action Wq 4 ^ for the onset of a plateau can be estimated 

to be bounded as 1.0 < W^f 1 < 1.2. As shown in the 
inset, these bounds are confirmed also by the behavior of 
the derivative dA(Co) / dCo, which is clearly larger than 
zero for Wo = 1.0, and becomes zero (within the limits 
of numerical accuracy) around Co — 0.5 for Wo = 1.2. 

For Co fixed the values of the prefactor A(Cq,Wq) 
as a function of Wq can be used to estimate via lin- 
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Co at the reservoir edge is, i.e., A(Cq, Wq > Ifj "') = 0. 
This value cannot be measured directly due to the unrea- 
sonably long simulation times needed to reach the asymp- 
totic regime in this range of Wo values, but the linear 
extrapolation of the available data, as shown in the inset 
of Fig.H», yields the estimate W (cot) ~ 2.3. 

In the range of low densities Co at the reservoir edge, 
the KMC results in Fig. ^indicate that there is a thresh- 
old value Cq TO ~ o.l for the density in the reser- 
voir edge below which, independent of the interaction 
strength Wo, there is no extraction of a monolayer: all 
the curves A(Cq) reach zero at a non-zero value of Co. As 
we will show below, this is a consequence of the condition 
(e) in the model, i.e., of the requirement that a move 
from r g r t toward a forward site r' is accepted only 
if there is at least another particle in the neighborhood 
\r' — r\ < r c . As we have mentioned before, this condi- 
tion is equivalent to requiring a minimum value Cl for 
the density on the advancing edge of the film. This den- 
sity Ci can be estimated as follows. For a fluid particle 
at (xq, ?/o ) G Ft (the filled circle in Fig.0 to move to the 
site (xo + 1, j/o ) (the open circle in Fig.EJ, in the shaded 
area of the disk of radius r c — 3 centered at [xq + 1, yo) 
there must be at least one more particle in addition to 
the one attempting the jump. The unshaded sector of 
the disk is due to the fact that by definition of T t the 
sites (x > xo,yo) are empty. Therefore, at least two of 



FIG. 4: (a) Dependence of the prefactor A(Co,Wq) on Co 
for several values of Wo. The inset shows the derivative 
dA(Co)/dCo as a function of Co in the cases Wo = 1.2 (O) 
and Wo = 1.0 (•); the dashed line in the inset corresponds to 
dA(Co) / dCo = 0. (b) Estimate of a "covering phase-diagram" 
in the Wo-Co plane showing the parameter range for which a 
macroscopic film is extracted from the reservoir or is not ex- 
tracted, respectively. The curve Wq CO "' shows the "covering 
- non-covering" separatrix described in the main text. The 
line is a guide for the eye, and the dotted line at Co < 0.2 is 
a heuristic extrapolation suggesting Wq C °"'(Co — > 0) ~ 1.4. 
The inset shows the dependence of A(Co = 1, Wo) on Wo. 
The line extrapolated to A(Co = 1, Wo) = is a linear fit of 
the last four points. In all four plots the symbols are KMC 
results; dashed and dotted lines connecting the symbols are 
guides to the eye. 



ear extrapolation, as shown in the inset in Fig. 01b) f° r 
the particular value Co = 1, the interaction strength 
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The data for W^° v) (C ) are 



(Co) at which for a given Co the prefactor van- 
ishes: A(C 0l W^ cov) ) = 0. 

shown in Fig. E{b). Since at a given Wo the prefac- 
tor attains its maximum for Co = 1, the dependence of 
^.(Co = 1, Wo) on Wo allows one to infer also an estimate 
of the interaction W^ cov ^ above which no macroscopic 

film is extracted from the reservoir whatever the density threshold value C, 



FIG. 5: Schematic drawing of the region around a point 
{xo, yo) (filled circle) belonging to the advancing edge Ft. The 
target destination, for the case of an advancing Ft, is denoted 
by an empty circle. The shaded area shows the domain in 
which at least one other site must be occupied so that the 
particle at (xo,yo) can move to (xq + 1, yo) in accordance 
with our KMC rule (e) (see Sect. II). Since (xo,yo) £ Tt, 
all sites (x > xo,yo) must be empty so that in the unshaded 
sector there are no occupied sites. 

the M = 25 sites in this shaded region are occupied, and 
thus Ci = 0.08 is an estimate for the minimum density 
on Tt ■ Since extraction of a film requires a moving edge 
T t , the density at the reservoir edge should be greater 
than Ci for spreading to be possible. This explains the 



(mm) 




0.1 observed in the simulations. 
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B. Asymptotic scaling 

We now turn to the analysis of the time dependence of 
the transversally averaged density profiles C(x, t) of the 
spreading monolayer. Since the time dependence of the 
advancing edge follows asymptotically X(t) ~ \ft in all 
the cases in which spreading occurs, it is natural to test 



if the density profiles C(x, t) actually scale as a function 
of the scaling variable A = x/^/Dot. In Fig. [5] we show 
density profiles C{x,t) for (a) W = 0.6, C = 1.0; (b) 
W = 1.4, C = 1.0; (c) W = 1.0, C = 1.0; and (d) 
Wo = 1.0, Co = 0.6 as functions of the scaling variable 
A = x/\/Dot, and as functions of x in the insets. 




X X 

FIG. 6: Density profiles C(x,t) for (a) Wo = 0.6, Co = 1.0; (b) W = 1.4, Co = 1.0; (c) Wo = 1.0, Co = 1.0; and (d) 
Wo = 1.0, Co = 0.6 as functions of the scaling variable A = x/\/Dot. The insets show these profiles as functions of x. The 
results correspond to t = ± T (O), ^ T (□), ^ T (O), and T (x) in (a), (c), and (d) and to t = T (O), 4 T (□), 7 T (O), 
and 10 T (x) in (b), where T = 2 x 10 6 . The space and time units are the same as those used in Figs. 121 II 



The results at late times show a very good data collapse 
and actually only the data corresponding to the earli- 
est time seems to show some significant deviations. This 
strongly suggests that in the asymptotic limit (t ^> 1, 
X(t) 3> 1) the density profiles can be described by a scal- 



ing function C(A = x/\/Dot; Wo, Co). Here we have ex- 
plicitly indicated the parametric dependence of the scal- 
ing function on the interaction strength Wo and on the 
density C at the edge of the reservoir (see Fig. |SJ). We 
note that for small Cq (like Cq = 0.6 in Fig. Eld)), one 
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observes deviations from scaling in the range of densities 
C(x,t) < 0.1. These deviations are most probably due 
to insufficient statistics, although they may also indicate 
that the true asymptotic regime has not yet been reached 
in the simulations. However, since even in this range 
there is a clear tendency of smaller changes in the shape 
of the density profile for increasing time t, it is reasonable 
to expect that the results in Fig. El are good approxima- 
tions for the corresponding scaling function C(A; Wo, Co). 

The scaled density profiles in Fig. El reveal three im- 
portant features. First, we have already noted that the 
change in shape of the function A(Cq) as Wq crosses 
1.0 < Wq < 1.2 signals a change in the spreading be- 
havior. As shown by the data in Fig. HJb), for large Wq 
the monolayer has an almost compact structure, and at 
the advancing edge there is a sharp transition from a 
large density to a small, almost zero, density. Therefore, 
in this range of interaction the spreading is accompanied 
by the emergence of a well defined interface between two 
phases. In contrast, at small Wo the density decreases 
smoothly from the value at the edge of the reservoir to 
zero and no jump in the density is visible. Thus W^p is 
a threshold value above which the attractive interaction 
is strong enough to support the build-up of an interface. 
In this sense, the change at W^ may be interpreted as 
the onset of a phase separation. 

The second point regards the shape of the density pro- 
file as a function of the parameter Co- For a weak attrac- 
tive interaction Wo (Fig. EI a )) or at small densities Co 
(Fig. Eld)), the density profiles resemble well the error 
function solution of a regular diffusion equation for non- 
interacting particles Q, Q • As we shall show in the next 
section, in this range such a description is not only qual- 
itatively but even quantitatively accurate. However, for 
larger values Wo (but still below W^) and for large Co 
(see Fig. EJc)), one finds the formation of a pronounced 
shoulder in the scaling function in the range of small A 
(i.e., A < 0.5), and thus a significant deviation from an 
error function solution. This shows that in this range of 
parameters the asymmetry in the jumping probabilities 
due to the attractive interaction between particles can- 
not be fully accounted for by an effective boundary force 
approach as in Ref . P, Q • Therefore one has to include 
explicitly this asymmetry into the description of the dy- 
namics in order to accurately capture the structure of the 
expanding film. Finally, interaction strengths above Wq 
have dramatic effects on the spreading behavior, leading 
to the emergence of interfaces (Fig.EJb)), and the simple 
description in terms of non-interacting particles breaks 
down completely. 

The third feature of the profiles to be discussed is the 
formation of a "foot" at the right end of the profile. The 
height of the foot is approximately equal to C\ and this 
is due to the fact that the density value on an advanc- 
ing edge T t cannot decrease below C\. The formation 
of the step implies that the fluctuations of the interface 
T t around the mean value X(t) are constant or increase 



in time slower than such that the width of the in- 
terface divided by \fDo~t vanishes in the long time limit. 
This sharp interface is occurring naturally due to the fact 
that the eventually large fluctuations are suppressed by 
blocking the advancing of isolated particles ahead of the 
film (see rule (e) in Sect. II), and thus the width of the 
interface would be expected to be of the order of the cut- 
off r c — 3 of the attractive potential and to be almost 
constant in time. Although the formation of the "foot" 
is a somewhat artificial feature introduced in the model 
by rule (e), the advantage of it is, as discussed in the 
Introduction, that it leads to the clear formation of an 
interface with its associated dynamics. 



V. CONTINUUM LIMIT 

A. Differential equation for the density and scaling 
behavior 

Neglecting all spatial and temporal correlations, i.e., 
assuming that averages of products of occupation num- 
bers rj(r;t) are equal to the corresponding products of 
averaged occupation numbers p(r;t) = (rj(r;t)), where 
(...) denotes the average with respect to the correspond- 
ing probability distribution V({r/(r;t)}) of a configura- 
tion {r](r] £)}, one can formulate the following mean-field 
master equation for the local occupational probability 
(density) p(r;t): 

M^ = -p(r;t) ^r'; t [l- P (r';t)} 

r f ,\r f — r\=l 

+ [l-p(r;t)] Uv>^v,tp{r';t) , (11) 

r', \r f — r\ — 1 

where 

U(r;t) = (U(r;t)) = -U £ S^ {12) 

r",0<\r"-r\<3 1 1 

is replacing U(r;t) in the definition for p(r — > r') 
(Eq. ©)■ 

As shown in detail in Appendix B, in the continuum 
space and time limit (At — > 0, a — ► 0, f2 _1 — ► 0, Do = 
f2a 2 /4 finite) of Eq. (|llfl . by taking Taylor expansions for 
p(r — > r') and p(r'\ t) around r and keeping terms up to 
second-order spatial derivatives of the density p(r; t) |37l 
EM l46| , one obtains the following nonlinear and nonlocal 
equation for p(r;t), 

d t p = D V[Vp + /3p(l- p)VU]+0(a 2 ) (13) 

Since the derivation of Eq. (|T3|> presented in Appendix 
B is not a rigorous proof (as we shall discuss below, such 
a proof appears to be extremely difficult to obtain) but 
rather a heuristic derivation in the spirit of Ref. |46| . 
several comments are in order before proceeding. The 
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only lattice gas system with long-ranged interactions for 
which it has been rigorously shown that Eq. (|13[) rep- 
resents the correct continuum limit at all temperatures 
is the hard-core lattice gas model with a Hamiltonian 
composed of short-ranged (on-site) repulsion and long- 
(infinite) ranged Kac potentials evolving via rates which 
satisfy detailed balance |47l I4H l49| . Recently, it has been 
argued that similar equations will also hold for systems 
with relatively short-rang ed interactions 50,51,521. The 
system discussed in Refs. [5(1 101 is a hard-core lat- 
tice gas model with attractive inter-particle interaction 
in the form of cither a finite range constant potential or 
of Morse potentials, a microscopic dynamics defined by 
nearest neighbor jumping rates depending on the energy 
at the departure site (Arrhenius dynamics) or on the dif- 
ference in energy between the two sites (Metropolis dy- 
namics), and a fixed density gradient imposed by fixed 
densities at the boundaries of the simulation box. As 
shown in Ref. [52J , the solutions of such continuum equa- 
tions compare fairly well with results of corresponding 
kinetic Monte Carlo simulations in cases where the inter- 
action range is greater than several lattice units. Typical 
values are 5-10 lattice units, depending on the dimension- 
ality of the problem, the type of potential, and the type 
of dynamics defined by the microscopic rates. Moreover, 
it has been shown by Leung |46| | that heuristic deriva- 
tions, based on formal Taylor expansions, of continuum 
equations from microscopic dynamics give very good re- 
sults both for the dynamics of Ising models with nearest 
neighbor interactions and Kawasaki rates and for that 
of driven lattice gases. Based on these results we assume 
that Eq. I|13|) is at least a good approximation for the the 
continuum limit of our system, although it was not rig- 
orously derived and although the range of the attractive 
potential in our problem is rather short (r c — 3). Finally, 
since the continuum limit for our system seems to be de- 
scribed by the same equation as that for the dynamics in 
systems evolving via rates preserving the detailed balance 
condition, one may conclude that indeed the deviations 
(calculated in Appendix B) from detailed balance in the 
rates defined by Eq. are small and do not carry over 
to the macroscopic scale. 



The constraint of a fixed density Co at the edge x = 
of the reservoir implies the boundary condition 



p(x = 0, y; t) = C . 



(14) 



As we have pointed out in Subsec. IV. A, the condition 
(e) in the model leads to a well defined interface and 
implies that for a spreading film the minimum density 
on the advancing edge is C\ > 0.08. In the absence of 
other additional constraints imposed by the formation of 
interfaces, i.e., for interactions Wq < \ and for large 
times, these considerations and the KMC results (see also 
Figs.UJa), (c), and (d)) strongly suggest that the density 
on the advancing edge X(t) can be considered as fixed 



and equal to Ci, leading to the boundary condition 

p(x = X(t),y;t) = C 1 . (15) 

In what follows, we shall use the value C\ = 0.11 obtained 
in the KMC simulations. We note in passing that the 
boundary condition Eq. (|15|) also naturally occurred in 
the "effective boundary force" model of Burlatsky et al. 
0,0, the expression of Ci in this case being C\ = 1 — /i, 
with p the ratio between the forward and the backward 
jumping rate for particles on the advancing edge. 

Since there are no boundaries along the y— direction 
and the boundary condition at the reservoir (Eq. I|14|) ') 
is independent of y, an important consequence of the y— 
independence of Eq. I|15|l is that the solution p(r;t) does 
not actually depend on y, which implies that the mono- 
layer is homogeneous along the y— direction, in agree- 
ment with the KMC results. Therefore one has to solve 
an effectively one-dimensional problem. The study of the 
occurrence of spontaneous transversal instabilities of the 
advancing edge would require to replace Eq. Ijlofl by a 
moving, transversally varying boundary condition. This 
might become relevant if the monolayer is driven by ex- 
ternal forces or encounters an obstacle. However, in view 
of the KMC results we have no reason to consider such 
effects for the present system. 

Although the reduced dimensionality is a significant 
simplification, Eq. i|13|) remains quite complex because 
it is nonlocal due to the term involving the interaction 
potential U(r;t). However, assuming that the density 
p{r; t) is a slowly varying function of the spatial coor- 
dinates (which certainly is a reasonable hypothesis ev- 
erywhere except near interfaces, see Figs. [5] and |SJ), the 
potential U(r;t) may be expanded as 



U(r:t) = -U 

r',0<|r'-r|<3 
r',0<\r'~r\<3 ' 



p(r':t) 



+ 0{a 2 ). (16) 



As discussed in Appendix B, the rotational symmetry 
of the lattice and of the factor \r' — r|~ 6 implies that 
the summation over r' will cancel the contributions of 
the first-order derivatives, and thus the leading gradient 
term does not appear in the expansion above. Because 
in the derivation of Eq. (|13fl only terms up to second- 
order spatial derivatives of the density have been kept, 
i.e., second-order in the lattice constant a, and a factor a 2 
has been already absorbed into the diffusion coefficient 
Dq, only the zeroth-order term in the expansion above 
will contribute. This leads to the local equation 

d t p = D V{[1 - g W Q p(l - p)] Vp} + 0(a 2 ), (17) 
where g = l r P 6 is a geometrical factor depen- 

l<\r\<r c 

dent on the lattice type (e.g., square, triangular, etc.) 
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and on the cut-off range of the potential. For the present 
case of a square lattice and a cut-off at r c = 3 one has 
g ~ 4.64. 

Rescaling the time as t — ► r = Z?o£ and defining an 
effective diffusion coefficient 



D e (p) = l-gW p(l-p), 



(18) 



Eq. (|17|l may be written in the usual form of a diffusion 
equation 



d T p = V[D e (p)V P }+0{a 2 ) 



(19) 



The functional form of D e (p) (Eq. I jlSjl l implies that 
for Wq > A/g there will be values pi of the density for 
which D e (pi) < (see Fig.CJ. For parameters such that 




itive argument allows an effective inclusion of correlations 
into the mean-field description and leads to a simple cor- 
rection to the mean-field value g = 4.64. The dynamics is 
possible only by jumps into empty locations. This means 
that the summation in g should include at most three 
contributions from nearest neighbor sites, giving g ~ 3.64 
and an estimate for the threshold interaction W^p ~ 1.1, 
in good agreement with the KMC results. Thus for the 
rest of the analysis we will use this corrected value of g. 

We now proceed with the analysis of the density pro- 
files for the asymptotic scaling limit. Since the solution 
p(r; t) depends only on x, Eq. i|19|) yields an equation for 
the transversally averaged density C(x, r) : 



d T C(x; t) = d x [D e {C)d x C(x; t)} + 0(a 2 ). 



(21) 



Introducing the scaling variable A = x/y/r leads to the 
following equation for the scaling solution (7(A) in the 
asymptotic limit t>l: 



XdC 
2~d\ 



d 

d\ 



+ O[(a/V^) 2 ]=0 (22) 



with the boundary conditions 

C(0) = Co 
C(A) = C x 



(23a) 
(23b) 



Since the solution of Eq. depends on whether Wo < 



W^ or Wo > W^ L ', we shall discuss these two cases 
separately. 



At) 



FIG. 7: Effective diffusion coefficient D e (p) (Eq. JTHJ) for 
gWo = 3 (upper curve) and gWo = 5 (lower curve), i.e., below 
and above the threshold value gW^ = 4, respectively. The 
values p~ and p+ indicate the range of densities pi for which 
Deifii) < 0, corresponding to instabilities in Eq. 119H . 

Wo < 4/g, Eq. is a proper diffusion equation, while 
for Wo > A/g instabilities are expected in the range of 
densities where D e (pi) < 0, i.e., for pi £ {p~ , p+) where 




It is known [47L l53| that these instabilities are discon- 
tinuities in the density profile ("shocks"), i.e., they cor- 
respond to the formation of interfaces, which is exactly 
what is observed in the KMC results. Thus the value 
for the threshold interaction strength Wq (introduced 
in Sec. IV) for which interfaces emerge is predicted by 
the continuum theory as Wq — i/g ~ 0.86. This value 
is significantly smaller than the lower bound estimate 
1.0 < Wq from KMC simulations, which is not unex- 
pected because of the mean-field character of the deriva- 
tion of the continuum equation. However, a simple, intu- 



B. Scaling solution for Wo < W { > 

For W < W (t) , in Eq. © the term C[(a/V7) 2 ] may 
be neglected and Eq. I|22|l together with the boundary 
conditions given in Eq. I|23|l is a well posed problem and 
thus admits a regular solution C{\]Wq,Cq). Although 
the solution cannot be found in closed form, the numeri- 
cal integration of Eq. lll'L'l) is straightforward. Results for 
small and intermediate values of the attractive coupling 
Wo and for several values of Co are presented in Fig. |H1 
For comparison, we also show results coresponding to the 
mean-field "effective boundary force" (EBF) approach 
subject to the same boundary conditions (Eq. (J2HI), for 
which the density profile is given by [54| 

a,„(A,= C „-( C „- Cl ,|S^, ,24) 

where erf(z) = f Q z dy e~ y is the error function. 

There is excellent agreement in all cases between the 
theoretical results from Eq. 11221) and the KMC results. 
Similar conclusions hold for all values of Co and Wq < 
1.0. (These results are not shown.) These findings offer 
additional strong support both to the assumption that 
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FIG. 8: Asymptotic scaling solution C(A) for (a) Wo = 
0.6, C(\ = 0) = Co = 1.0, 0.8, 0.4 and (b) W = 1.0, C(X = 
0) = Co = 1.0, 0.8, 0.4 with A = x/y/Dvi. Shown are 
theoretical mean-field results from Eq. 12211 (solid lines), re- 
sults of the "effective boundary force" (EBF) theory from 
Eq. 12411 (dashed lines), and corresponding KMC results at 
time T = 2 x 10 6 (O) (assumed to be close to the asymptotic 
limit). 



Eq. IL'l'l) is an accurate description of the continuum limit 
for W < W (t) and to our heuristic correction g ~ 3.64 
for taking correlations into account (see the paragraph 
preceding Eq. iJHJ). When compared to the EBF results 
from Eq. i|24|). we see that, as expected, at low densi- 
ties of the reservoir (see, e.g., the curves in Fig. |S] cor- 
responding to Cq = 0.4) the predictions of Eq. 11' I'll and 
the EBF results are almost identical because the mono- 
layer is dilute and thus the particle-particle interactions 
are less effective. For similar reasons, at low values for 
the strength of the inter-particle attraction or at high 
temperature, i.e., when Wo is small, the EBF description 
performs well even for high densities (see in Fig. Efb) 
the curve corresponding to Co = 0.8). However, at very 
high densities of the reservoir or for large values of the 
attractive coupling Wo, there are significant discrepan- 
cies even in the qualitative behavior between the EBF 



predictions and the simulations results. In particular, 
the formation of a "shoulder" in the case in which Wo is 
large (see in Fig.UJa) the curves corresponding to Co = 1 
and Co = 0.8) is remarkably well reproduced by the theo- 
retical curve obtained from Eq. (|22|l . but it is completely 
missed by the EBF solution (Eq. J33J). Therefore we 
conclude that even in this case, i.e., below the threshold 
value Wq for interface formation, the inter-particle at- 
traction has to be explicitly included into the model in 
order to obtain a correct prediction for the mass distri- 
bution inside the monolayer which is extracted. 

The above results should also be discussed in the con- 
text of the similar work in Refs. [5l| and [S^] mentioned 
in the beginning of this subsection. We have emphasized 
that the derivation of the continuum limit is mean-ficld- 
like in character, and that only after correcting for corre- 
lations, i.e., after adopting the improved value g ~ 3.64, 
the continuum limit accurately predicts both the thresh- 
old value Wq^ for the interaction coupling and the scaled 
asymptotic density profiles C(A;Wo,Co). Such a cor- 
rection has not been included in the similar continuum 
equations discussed in Refs. (5l| and [5^], and we suggest 
that this explains the discrepancies observed by the au- 
thors in the case in which the range of the inter-particle 
potential is short (see, e.g., the density profiles corre- 
sponding to cut-off ranges r c = 2 and r c — 5 in Fig. 4(a) 
in Ref. [52] )• For the longer-ranged potentials (r c > 5) 
used in Refs. [HJ and further than nearest neighbors 
contribute significantly to g, and thus the exclusion of a 
nearest- neighbor term becomes relatively less important, 
which explains the good agreement obtained in the cases 
r c > 5 without any correction included. 

Before proceeding to the case Wo > Wq \ we would 
also like to briefly comment on the connection between 
our above results, the experimental results for the pre- 
cursing films of Pb on Cu(lll) reported in Ref. jE^, and 
the MD results for precursing films of Ag on Ni(100) 
presented in Ref. |56j| . The density profiles measured ex- 
perimentally (see Fig. 3(b) in Ref. |5^) and in the MD 
simulations (see Fig. 3 in Ref. |5(j) for the 2D spreading 
of Pb or Ag films show a striking resemblance with the 
ones we have obtained in the KMC simulations. More- 
over, by assuming a macroscopic diffusive dynamics de- 
scribed by an equation of the same form as the one de- 
rived in Eq. (|21|l . effective diffusion coefficients D e (C) 
have been obtained from the density profiles, and the 
data shown in Fig. 3(a) in Ref. |55j are in qualitative 
agreement with a quadratic dependence for D e (C) as in 
Eq. (|18|l . Therefore, the very simple microscopic model 
we discussed seems to capture the essential features of the 
dynamics in these cases. Moreover, this suggest that this 
form of D e (C) for Pb on Cu(lll) is not necessarily due 
to surface alloying - a mechanism which is not included 
into the dynamics of our model - but is rather already a 
consequence of the interplay between the concentration 
gradients and the inter-particle interaction, which leads 
to "jamming" and thus significantly slows down the dif- 
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fusion for large values of the density C(x,t). 



conditions expressed in terms of n(C): 



C. Scaling solution for Wo > W^ t] 

We now turn to the discussion of Eq. ill'L'l) for the case 
Wo > Wq . Because in this case the effective diffu- 
sion coefficient D e (C) becomes negative within a range 
of densities, the problem is known to be mathemati- 
cally ill-posed and to lead to discontinuities ("shocks") 
in the long-time limit if the small terms 0[(a/^/i) 2 ] (see 
Eq. (|22fl ) are set to zero [H^- For this problem the exis- 
tence and uniqueness of a "weak" solution C(X) ("weak" 
in the sense that C(X) has a discontinuity at a point 
A = A s but satisfies Eq. (12211 for A ^ A s ) have been re- 
cently addressed by Witelski jS^, EH using singular per- 
turbation methods. We will use here directly the explicit 
construction of the shock solution derived in Ref. (53 for 
the case in which the term 0(a 2 /y/t) is proportional to 
d^.C(x,t), the details of the calculation being presented 
in Appendix C. 

Defining 



(25) 



H{C)= / dC'D e {C) 
Jo 

Eq. H21|) may be rewritten as 

d T C = d 2 x p{C) + 0{a 2 ) , 



(26) 



i.e., it has the form of a diffusion equation for C(x,t) 
with a mobility M = 1 and a "chemical potential" fi(C). 
Moreover, as we will discuss below, the values of fJ.(C) 
across the discontinuity satisfy conditions which are sim- 
ilar to those determining the equilibrium liquid-vapor co- 
existence line in the van der Waals-Maxwell mean-field 
theory of liquid-vapor transitions. Because of these sim- 
ilarities, in what follows we shall informally denote /i(C) 
as chemical potential. 

Following Refs. [S3, E3, we look for a weak solution 
of Eq. subject to the boundary conditions given in 
Eq. 123(1 . in the form of a "shock" defined as 



(7(A) 



Ct{\), 

a (A), 



A < A s 
A > A s 



(27) 



where Ci(X) and C r (X) satisfy Eq. 11221) in the intervals 
[0, A s ) and (X S ,A = X(t)/y/i], respectively, subject to 
the boundary conditions 



Ce (0) = Co, Cg(X s ) — Cm, 
C r (X s ) = C m < Cm, C r (A) = C\, 



(28) 



respectively. As discussed in Appendix C, the singular 
perturbation analysis of Eq. (12211 implies that the values 
Cm and C m of the density at the left and the right of the 
shock, respectively, are determined from the following 



Cf, 



H(C M ) = MCm), (29a) 
dC [m(C") - fi(C M )} = 0, (29b) 



i.e., continuity of the "chemical potential" n(C) across 
the shock and a Maxwell equal area-rule across the shock, 
as mentioned at the beginning of this subsection. Solving 
Eq. (|29p. we find that the only solution satisfying the 
condition Cm > C m is 



Cm 
C m 



1 + ll i-J- 

2 2 V gw 



V3 
2 



1 - 



gW 



(30a) 
(30b) 



Comparison with the similar expressions for C^ , where 
D e (C±) = (Eq. GDI), shows that for all W > W (t) one 
has C m < C~ < C+ < Cm, thus the shock occurs both 
above and below the interval corresponding to unstable 
states. For the states corresponding to densities C G € = 
(C m , C~)U (C+ , Cm) the effective diffusion coefficient is 
positive, but the density gradients are very large in the 
long-time limit and the state becomes part of the shock; 
thus densities C £ €. correspond to metastable states. 

The last unknown, the position A s of the shock, is ob- 
tained from the conservation of mass. In an infinitesimal 
time interval St, the displacement 5s = x s (t+St) — x s (t) 
of the position x s — X s y/r of the shock leads to an 
increase (Cm — C m )Ss in the mass inside the stripe 
x s (r + St) — x s (r) . This should be equal to the net mass 
transfer ST[j(x s ) — j(x s + 6s)], where the mass current 
j(x) — —d x fi(C) (see Eq. lHHJ) is discontinuous at x s . 



1/2 dC 



Since 5s /St = (1/2) X s t-^ 2 and d x C = t 
obtains the following expression for the position A s of 
the shock 



D e (C 



M 



A s 



dC 
~dX 



D e (C n 



dC 
~dX 



Cm — C„ 



(31) 



We note that the above result can be also obtained via a 
direct integration of Eq. (|22|l across the shock, i.e., from 
A = A s — £ to A = A s + £ in the limit £ — * 0, using 
for the density profile there the approximation by a step 
function C(A) = Cm ~ (Cm ~C m )H(X~ X s ), where H(x) 
is the Heaviside function (H(x < 0) = 0, H(x > 0) = 1). 
It is important to note here that for sufficiently large 
values Wo of the attractive interaction the density C m 
may become smaller than C\. Since the density at the 
advancing edge cannot be smaller than C\, in this case 
the branch C r (X) disappears and the shock position is 
obtained by setting C m = in Eq. I|31|) . 

Once C m and Cm are known, Eqs. ffify. iffify . and J3U 

can be, in principle, solved for the corresponding quanti- 
ties Ci(X), C r (X), and the position X s of the shock. Since 
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Eq. (1221 cannot be solved in closed form, the above sys- 
tem of equations has to be solved numerically. Such a 
numerical solution is shown in Fig. for the cases (a) 
W = 1.2, C = 1.0, for which C m ~ 0.25 > C x , and (b) 
Wo = 1.4 and C = 1.0, for which C m ~ 0.098 < Ci. 
It can be seen that the agreement between the theoreti- 
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FIG. 9: Asymptotic scaling solution C(A) for (a) Wo = 1.2, 
Co = 1.0, and (b) Wo = 1.4, Co = 1.0. Shown are theo- 
retical results obtained from Eqs. 1221 . (12811 . and 1311 (solid 
lines), and corresponding KMC results at time T = 2 x 10 7 
(O) (assumed to be close to the asymptotic limit). The dot- 
ted line is a guide to the eye. The dashed lines indicate the 
corresponding values C m and Cm from Eq. 130fl . C& from 
Eq. 1201 where D e (C^) = 0, which determines the onset of 
the density range leading to instabilities (see Fig.0, Ci from 
the boundary condition Eq. 123bl . and the position X s of the 
discontinuity given by Eq. 1311 . 

cal asymptotic "shock" solution and the KMC measured 
density profiles is good for the large value Wo = 1.4, but 
it is not so good in the case Wo = 1.2. This is very likely 
due to the fact that in the latter case the simulation has 
not yet reached the true asymptotic regime, while for 
Wo = 1.4 the approach to the asymptotic shape is faster 
because the low density branch C r is suppressed. In both 
cases the KMC results confirm the value Cm as the onset 
of large density gradients, and there is good agreement 
between the theoretical prediction and the simulations 



in the range of densities C > Cm- This also supports 
the above conclusion that the discrepancies in the range 
C < Cm are due to simulations times which are not large 
enough. 

The jump Cm — C m in the density at X s explains the 
formation of the plateau (for Wo > W^) in the depen- 
dence of A(Co, Wo) on Co- if the density Co at the reser- 
voir is within the range C m < Co < Cm, in the imme- 
diate vicinity of the reservoir edge the density drops to 
C m and in the long-time limit the extraction of the film 
proceeds effectively as if the reservoir density would have 
been C m . Also, since 1/2 — C m = — 1/2 + Cm, it fol- 
lows that the plateau should be symmetric with respect 
to C = 0.5; indeed the KMC data in Fig. Ufa) exhibit 
this symmetry (as long as Wo is such that C m > 0.1). 
Moreover, since the density must satisfy C < 1, one 
may conclude that for interaction values Wo such that 
Cm > 1 the extraction of a monolayer is no longer pos- 
sible. This implies that the exact value for the upper 
limit of the interaction wjf ^ above which no macro- 
scopic film is extracted from the reservoir is given by 
6/g ~ 1.65. This value is significantly below 



W, 



(cov) 



the value W, 



(cov) 







W, 



(cov) 







(Co = 1) — 2.3 extracted from 
the linear extrapolation of the KMC data (see Fig.^Jb)), 
the discrepancy very likely reflecting that the KMC simu- 
lations have not yet reached the true asymptotic limit (or 
that a linear extrapolation is not appropriate) . Thus it is 
to be expected that in the range Co > 0.85 the separatrix 
^y(cov) g i 10wn j n Fig.^Jb) significantly overestimates the 
correct curve. 



VI. SUMMARY AND CONCLUSIONS 

Using kinetic Monte Carlo (KMC) simulations and a 
non-linear diffusion equation within the continuum limit, 
we have studied a lattice-gas model with interacting par- 
ticles for the two-dimensional spreading on homogeneous 
substrates of a fluid monolayer which is extracted from a 
reservoir (Fig. We have obtained the following main 
results: 

1. The two-dimensional KMC simulations confirm the 
time dependence X(t — » oo) = A\ft of the spreading, 
where X{t) is the average position of the advancing edge 
of the monolayer at time t, and reveal a non-trivial de- 
pendence of the prefactor A on the strength Uq of inter- 
particle attraction and on the fluid density Co at the 
reservoir (see Figs. [21 01 and QJ. A careful analysis 
of this behavior has allowed us to identify, in terms of 
Wo = Uo/ksT, a transition point W ~ 1.1 associ- 
ated with the occurrence of interfaces inside the extracted 
monolayer, and to estimate a "covering phase-diagram" 
in the Wo~Co plane (Fig. 0} together with a "covering 
- non-covering" separatrix W below which a macro- 
scopic film is extracted from the reservoir, while above 
y^(cov) ^ ^ extracted. 
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2. The asymptotic (i.e., at long time and large spa- 
tial scales) transversally averaged density profiles C{x, t) 
measured in the KMC simulations exhibit a scaling be- 
havior as function of the scaling variable A = x/^/D$t, 
where Dq is the one-particle diffusion coefficient on the 
bare substrate (Fig. EJ) . They clearly show that for this 
model the density in the extracted monolayer is not spa- 
tially constant, in contrast to the predictions of other 
theoretical models mentioned in the Introduction. This 
provides an unambiguous - and otherwise difficult - way 
to experimentally discriminate between the various theo- 
retical models proposed. Moreover, the simulations show 
that the present model predicts qualitatively different 
structures for the experimentally accessible density pro- 
files above and below the threshold value W (see Figs. [3 
and[SJ|, in particular the formation of sharp interfaces in- 
side the extracted monolayer for Wq > 

3. The asymptotic, scaled density profiles C(X) have 
been analyzed within a continuum limit with the cor- 
responding non-linear diffusion equation derived from 
the microscopic master equation. Within this approach 
we have included the effect of correlations in an effec- 
tive manner into the standard mean-field description by 
adapting the value of the integrated attractive interaction 
to account for the presence of empty nearest-neighbor 
sites (see Fig. |5J ■ This leads to an excellent agreement 
between the theoretical predictions based on the contin- 
uum limit and the KMC results both for the value Wq^ 
and for the scaled density profiles (Fig. [HJ) . Additionally 
we have shown that, even below the threshold value Wq^ 
for interface formation, the inter-particle attraction has 
to be explicitly included into the model in order to obtain 
correct predictions for the mass distribution inside the ex- 
tracted monolayer. The formation of the interfaces in the 
range Wq > Wq has been related to instabilities of the 
diffusion equation associated with densities for which the 
corresponding effective diffusion coefficient becomes neg- 
ative (Fig. EJ. We have constructed the corresponding 
discontinuous density profiles ("shocks") and critically 
compared them with the KMC measured ones (Fig. |UJl. 
Based on the results of a singular perturbation analy- 
sis, we have obtained a good estimate Wq° ~ 1.65 for 
the upper limit of the interaction above which no macro- 
scopic film is extracted from the reservoir. 

Finally, we comment on the connection between this 
model and experimental systems. As briefly discussed 
in Subsect.V.B, the present model appears to provide a 
successful description for the diffusion of solid metals on 
metal surfaces as studied in Refs. |55| and |56j. We have 
found a qualitative agreement between the experimental 
results in the case of diffusion of Pb on Cu(lll) [55| and 
our theoretically derived density profiles and effective dif- 
fusion coefficient. It seems to be promising to explore 
quantitatively the applicability of the present model for 
such metal on metal systems. To this end the experimen- 
tal setup described in Refs. |5!| and j5|| would have to 
be modified in order to have straight instead of circular 



spreading geometries and a deposit-substrate combina- 
tion chosen such as to avoid surface alloying effects. 

As noted in the Introduction the experiments with flu- 
ids performed so far deal with polymer oils. As long as 
the entanglement of the polymer chains is not important, 
one may consider a coarse-grained description in which 
the chain is replaced by an effective particle of the size of 
the corresponding radius of gyration and only the motion 
of the center of mass is considered. Although the motion 
of these effective particles might not resemble simple, ac- 
tivated jumping processes so that the microscopic model 
description is not directly applicable, it is reasonable to 
expect that the macroscopic evolution will be diffusive. 
Therefore, one may expect that the continuum limit of 
the present model can be used to describe the spread- 
ing behavior observed in experiments with polymer oils, 
but the macroscopic parameters entering into the diffu- 
sion equation should be regarded as fit parameters and 
not quantities calculated from microscopic dynamics as 
considered here. 

In order to obtain a direct, quantitative test of the 
present theoretical predictions for precursor liquid films, 
new experiments would have to be performed using sim- 
ple liquids chosen such that they have a spreading rate 
large compared with the evaporation rate. This should 
be combined with observation techniques chosen such 
that the density profiles, and not only the spreading rate, 
could be measured, which would require an in plane (lat- 
eral) resolution in the order of few nm for the case of sim- 
ple liquids, i.e., several lattice constants, and in the order 
of 10-50 nm for the case of polymer oils, i.e., several inter- 
"effective" particles distances, because the density varia- 
tions are expected to occur on larger length scales. One 
technique which possibly may fulfill these requirements 
is reflection interference contrast microscopy |59j , assum- 
ing that the microscope objective may scan the area of 
interest (of the order of mm 2 ) in times sufficiently small 
compared to those on which the density profile changes. 
The technique has been used before in studies of (equilib- 
rium) wetting properties on micropatterned solid surfaces 
|59| , and already at the time of its first implementation 
a lateral resolution of at least 200 nm (see Fig. 13(b) in 
Ref. ) combined with a normal resolution of the order 
of 1 nm has been achieved. 
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APPENDIX A: KINETIC MONTE CARLO 
METHOD 



APPENDIX B: HEURISTIC DERIVATION OF 
THE CONTINUUM LIMIT 



In this section we discuss in some detail the variable- 
ste p co ntinuous-time kinetic Monte Carlo algorithm |4ll 
l43| that we have used. The main idea is to consider 
the sequence of independent, uncorrelated events repre- 
sented by jumps of particles away from the wells in which 
they were residing. Each of these events has an identical 
time- and environment-independent rate f2 as shown by 
Eq. J5J, in contrast to the location-dependent rates of 
particular transitions r — * r' (Eq. J3J). 

Consider the system at time t when there are N par- 
ticles in the film (a; > 0) and an event just occurred. 
Since the attempts of any particle to leave its well are un- 
correlated to similar events of other particles, and since 
for each particle the rate for a successful jump is Q, for 
each one of the particles the probability that until time 
t' > t no successful attempt occurs is P\{t) — exp(— Or), 
where r =t' — t. Therefore, since the jumps are uncorre- 
lated, the probability that none of the N particles expe- 
rienced a successful attempt in r is Pn{t) = [Pi(t)] n = 
cxp(-NHt) and the probability that the first successful 
jump will take place at if will be given by 



P = NVLP n (t) = iV^exp(-iV^T). 



(Al) 



Thus the time interval r between successful jumping 
attempts (between "events") is a random variable dis- 
tributed according to Eq. (|A1I) . Since all the N particles 
have identical rates f2 for events, the probability for a 
certain particle to be the one undergoing the jump is 
fi/(iVf2) = 1/iV, i.e., the particle to jump is selected at 
random. Let us assume the selected particle is at lo- 
cation r. There are z = 4 nearest- neighbor locations, 
and thus z — 4 possible realizations of the jump; the 
one to be actually attempted is selected according to 
the probability defined by Eq. Q. Specifically, calling 
the four probabilities pi, . . . ,pi, with pi corresponding 
to the jump (x,y) — > (x + l,y) and the others being 
indexed counter clockwise, one compares the successive 

k=j 

sums so = 0, Sj = ^""^ pfc, j = 1,2,3,4, with a random 

fc=i 

number v £ [0, 1] and selects pk for which Sk-i < v < Sk- 
As described in the text, the jump takes place if the se- 
lected destination site is empty, and is rejected if the 
destination site is occupied. 

We note here that, as shown in Ref. |43|. incrementing 
the time between events using intervals generated accord- 
ing to Eq. IjAlll and not a constant time interval equal to 
the average time l/(NSl) between events, like in a classi- 
cal Monte Carlo simulation, is essential in assuring that 
the simulated time is the correct real time, and thus that 
the simulations capture the correct time development of 
spreading. 



The (mean-field) master equation for the local occupa- 
tional probability (density) p(r;t) is given by 



dp{r;t) 
dt 



= -p(r;t) uj r ^r'-,t[l - p(r';t)} 

r' , If' — r| — 1 

+[l-p(r;t)] Y u r '~,r;tp(r';t), 

T"',|t*' — r\ — 1 



(Bl) 



where 



LU r — >r '-t — Q 



expj 




;t)-U(r' 


;*)]} 


E 

f f 1 7 ' ' f 1 = 1 


J/3 
cxpj- 


[U(r;t)- 


tf(r';t)]} 



and 



U(r;t) = -U Y, 

r", 0<|r'-r|<3 



(B2) 

P(r";t) 



We consider a two-dimensional regular lattice of coor- 
dination number z and lattice constant a and choose the 
orthogonal x — y coordinate system such that the a;— axis 
is along one of the lattice directions. For a given site r 
we index the nearest neighbors as r', j = 0, 1, 2, ...z — 1, 
where j = is chosen such that r' — r is parallel to the 
x-axis and j runs in counterclockwise direction. Denot- 
ing the angle formed by the vector r'- — r with the cc-axis 
j 

as <f)j so that </>,■ = 2tt— , the components x\ and y' of 

z J J 

r'j — r are given by x'j — acos(4>j) and y'j = asin(<3 
The following relations are satisfied by the angles <j)j 
and will prove to be useful for the rest of the calculation: 



2-1 



5>n fc (^) = 0, ^cos fc (0,)=O, 
k odd, < k < z, 



2-1 



z 
2' 



^sin 2 (0j) - -, J]cos 2 (^ 

3=0 3=0 

2-1 2-1 

^sin(2^)=0, 5^cos(2^)=0, 

3=0 3=0 



(B4a) 
(B4b) 

(B4c) 



Defining 



Sh{r,r' j ;t) = h{r' j ;t)-h(r;t), 



(B5) 



where h(r;i) is any of the functions p(r;t), U(r;t), or 
products of them, expanding h(r'y,t) near r, and sum- 
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ming 5h(r, r'-; t) over r'j one obtains 



3=0 



9h ... dh z r - 1 . 

— ^cos(^) + — ^sm(^) 



a 

IF 



^ £ cos2 (^) + 2 g^gZ 22 sm(^ ) cos(^ ) 
3=0 ^ 3=0 



<9 2 /i z ~ 

■— ^Tsin 2 ^) 



(B6) 



In the relation above, h = /i(r; t) and the derivatives are 
evaluated at r. Replacing the corresponding sums by the 
results in Eq. (|B4|> . it follows that 



a , and therefore the other terms on the RHS of Eq. (|B10|) 
will lead to contributions proportional to a 4 and higher 
orders. Therefore, up to contributions which are of sec- 
ond order in the lattice constant in the master equation, 
Eq. IjBlOfl may be rewritten as 



n 

uj r ^ r '-t — — exp 

z 



(Bll) 



This means that the deviations of the rates oj r _» r ' ; t from 
detailed balance, which according to Eq. (|B10|) are of 
second order and higher in the lattice constant, in the 
equation corresponding to the continuum limit contribute 
with terms of fourth order and higher in the lattice con- 
stant. These terms become negligible in the limit a — > 0. 



z-l 



T£Sh(r,r' j ;t) = -?-V 2 h + 0(a 4 

3=0 



(B7) 



Straightforward algebra allows one to derive from the 
definition l|B5|) the following additional useful relations 
involving a second function f(r;t): 

f(r' j ;t)Sh(r,r' j ;t) = 5(fh)(r,r' j ;t) 

-h(r f ,t)5f(r,r' r ,t) (B8a) 
Sf(r,r>;t)5h(r,r' r ,t)=8(fh)(r,r>;t) 
-M*j;<)*/(r,r5;t) - /(r i; t)fl»(r,r<;t). (B8b) 

Assuming that U(r; t) varies slowly on the scale of the 
lattice constant so that (3SU(r,r';t) <C 1, one has the 
expansion 



exp 



■^SU(r,r' k ;t) 



exp 

3=0 



-%SU(r,r' j; t) 



exp 



-Z6U(ry k ;t) 



{l - P SU(r, rj; t)/2 + [(3 SU(r, rj; t)/2] 1 



3=0 



Thus using Eqs. (B^ and p38b|) one obtains 



•■} 

(B9) 



The expression (IB11|) may be now expanded in terms 
of powers of (36U(r,r';t): 



where 



;t -P0+P1+P2 + 



Po = — , 

z 

Pi = -^SU(r,r>;t), 
2z 



P2 



n/3 2 



[SU(r,r';t)f 



(B12) 

(B13a) 
(B13b) 
(B13c) 



Note that the expression i|Bll|) implies 



ClV— >TVt 



exp 
exp 



P 



5U(r',r;t) 



+^su( r y,t) 



(B14) 



and thus its expansion in powers of (3SU(r, r'\ t) is 

kV->r;t ~ Po - Pi + P2 + ■ ■ ■ ■ (B15) 



We now will compute separately the contributions of 
these terms to the RHS of Eq. I|B1|I . For the contribution 
due to po one has 



■ exp 



P 



6U(ry k ;t) 



sil + t-f [V 2 U - P(VU) 2 ] +0(a 4 ) 



(BIO) 



where as before U = U(r;t) and the spatial deriva- 
tives are evaluated at r. As we shall show below, the 
zeroth-order term in the above expansion (Eq. (IB10|) ) 
contributes to the master equation already in the order 



z-l 



RHS ( o) = -X>(r<,i)[l-p(r,t)] 

Z 3=0 

- p(r,t)[l~p(r' v t)}} 



3=0 

{la 2 



-v 2 P + no(a 4 ). 



(B16) 
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For pi one has 



rhs (d = ^E^^H^^ 1 -^'*)] 

Z 3=0 

+ p(r,t)[l-p(r' v t)}} 

3=0 

z-1 
3=0 



2z 
2z 



(B17) 



Using Eqs. (|B7|I and (|B8a|l to replace the two sums in 
the expression above, one obtains 

pn a 2 



RHS (1) = {(1 - 2p)[V 2 (pU) - UV 2 p] + P V 2 U} 

QO(a 4 ) = V[p(l - p)VU] + flO(a 4 ) . (B18) 



Finally, for p2 one has 



RHS 



(2) 



4z 

3=0 

p(r,t)[l-p(rj,*)]} 

/3 2 Q ' ' 
4^ 



^[<5{/(r,^.;t)] 2 {p(^,t)[l-p(r,0] 



^[SUiry^fSpiry^t). (B19) 



Using repeatedly Eqs. (|B8a|) and (|B7() in the above sum, 
one obtains 



RHS 



(2) 



(3 2 Qa 2 
4 



[V 2 (pU 2 ) - pV 2 U 2 - U 2 V 2 p 



- 2UW 2 (pU) + 2UpW 2 U + 2U 2 V 2 p] 

+ nO{a A ) = ilO(a 4 ) . (B20) 

It is easy to see that higher order terms, p^,..., will 
contribute with terms which are at least of the order a 4 , 
and thus the expressions l|B16() and (|B18|) are the only 
terms relevant for Eq. i|Bl|) . 



Collecting the terms and passing to the continuum 
limit a — ► 0, SI -1 — > such that Do = £la 2 /4 stays 
finite, one arrives at the result given in Eq. in the 



d t p = D V{Vp + 13 [p(l - p)VU]} + 0(a 2 ) . (B21) 



APPENDIX C: DERIVATION OF THE SHOCK 
SOLUTION 

Following Ref. we start from equation Eq. lt27)|) 
written in terms of the scaling variable A as 



l.dC d 2 
2 ~d\ = dA 2 



Q 



d 4 C 

Ix 4 ' 



d 2 C 

d\ 2 



(CI) 

where the function Q is a linear combination of fourth- 



order derivatives terms of the form ^jj, (w&j , 

d 2 (UC) ~ " 

d\ 2 



... The region of interest, A € [0, A], nat- 
urally decomposes into the region near the interface, 
A s — h(e) < A < A s + h(e), and the outer region 
| A — A s | > h(e) with e = a/\fr; h(e) is a smooth function 
such that h(e — > 0) = (which ensures that in the long- 
time limit the width of the interface becomes negligible) 
and lim c ^o — * oo, i.e., it is assumed that the decrease 
of the width is slower than e. In the outer region, the so- 
lution Ci_ r (X) is a slowly varying, smooth function of A, 
and the terms proportional to e 2 in Eq. IjClfl are negli- 
gible. In contrast, in the inner region the gradients are 
very large, and the fourth-order terms become relevant. 

In order to obtain the shock structure, we change to 
the "stretching" variable £ = (A— A s )/efor |A— A s | < h(e) 
and look for a smooth, strictly decreasing solution C(£). 
In terms of £ Eq. i|2ti|) turns into 



which in the limit e 
e- approximation : 



2 Q 



d 4 C 



d 2 C 

~dC 



(C2) 



leads to the zeroth order in the 



Q 



d 4 C 



d 2 C 

d( 2 



(C3) 



Since in the limit e — > the inner region —h(e)/e < ( < 
h(e)/e is mapped into — oo < C < +°° an d at the ends of 
the shock region the inner solution must match the outer 
solution Ci tr (\), Eq. (|C3(I should be solved subject to the 
boundary conditions 

lim C(() = C M , km C(0 = C m , (C4) 

C — > — oo C — >+oo 

which implies also that all the derivatives (£) of the 
smooth inner solution C(£) tend to zero as C — > ±oo. 

Although the function Q may be computed explicitly 
by using Eq. I|B7|) , 



z-1 

E 

3=0 



5h{ry-t) 



za 



V 2 h 



3za 



V 4 h + 0(a 6 ), (C5) 
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and following similar steps as in calculating the terms 
proportional to a 2 in Appendix B, the result is very com- 
plicated and a singular perturbation analysis appears to 
be extremely difficult, if at all possible. However, one 
may argue that the term d 4 C is always relevant in the 
region of the shock because it is associated with an inter- 
face contribution Ti = J dx[WC(x, t)] 2 to the free energy 
functional T — T% + ■ ■ • of a dynamical Cahn-Hillard the- 
ory of phase separation, d t C = V [M(C)V ($§)] (where 
M(C) denotes the mobility) [5^ . Therefore, all the other 
terms that are relevant should be of the same order as 
d*C. This leads us to the approximation 



Q 



d 4 C 

d( 4 ' 



d 2 C 

d( 2 



(C6) 



with q a constant or a very slowly varying function of A, 
e -g-> q(0 — q(^s +eC)- I n this case, Eq. I|C3|> reduces to 



d?_ 

d( 2 



d 2 C 



o. 



A first integration leads to ^ /^(C) + q 



because both ^//(C) and ^§ are zero at £ 
and thus a second integration yields 



M (C) 



d 2 C , 



where b is an integration constant. Since limf_>± 



(C7) 
const = 

— > ±00, 

(C8) 



d-'C 



0, one obtains /i[C(£ — > ±oo)] = 6, i.e., the requirement 
of continuity of //(C) across the shock (Eq. (|29a|l in the 
main text): 



m(Cm) = m(C,„) ■ 
Since ^ ^ (except at infinity) and b 



Eq. I|C8|I may be rewritten in the form 

2" 



q d 
2dC 



dC 



\m(c m ) - M (C)] 



(C9) 
(CIO) 



which leads to 



dC[/x(C M )-M(C)] 



00 <f 



dC 



(Cll) 

i.e., the equal area-rule for (x(C) (Eq. (|29bjl in the main 
text). 



Finally, we remark that all the details of the calcula- 
tion, as well as the main results (Eqs. I|U9|I and (|Cll|l ). re- 
main unchanged if the corrections Q would have the form 

Q[cM,(ceo) 2 ,...] = ^p[c,c( 2 ),(c«) 2 ,...], withp 

a linear combination of terms of second order spatial 
derivatives satisfying lim^-i-ooP — ► [6T| . and thus it 
seems reasonable to assume that in general it is only the 
inner structure of the shock that depends on the partic- 
ular form of g mm. 
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